Estudio de lugares de observación óptima para el asteroide 99942 Apophis

Agustín Vallejo Villegas


Artículo #1

Astrodinámica

Universidad de Antioquia

2020-1

Librerías

In [1]:
"LIBRERÍAS BÁSICAS"
import numpy as np
import pylab as plt
import pandas as pd
import matplotlib.animation as animation
from matplotlib import rc
rc('animation', html='jshtml')
In [2]:
"SPICEYPY"
!pip install spiceypy
import spiceypy as spy
Collecting spiceypy
  Downloading https://files.pythonhosted.org/packages/fd/48/ad651637d85e7a3e517025babb1cb478cc20a18a4a6ad1a0f1463a8a5672/spiceypy-4.0.0.tar.gz (265kB)
     |████████████████████████████████| 266kB 4.2MB/s 
Requirement already satisfied: numpy>=1.17.0 in /usr/local/lib/python3.6/dist-packages (from spiceypy) (1.19.5)
Building wheels for collected packages: spiceypy
  Building wheel for spiceypy (setup.py) ... done
  Created wheel for spiceypy: filename=spiceypy-4.0.0-cp36-none-any.whl size=1830970 sha256=a443c1eef86f79f180f12773452dd77926f9585551643d56eea14e80e756ff1d
  Stored in directory: /root/.cache/pip/wheels/3d/ec/36/9d76db48ce8f8e50de3875c3e6650aefb2f9c0888dd2dca15c
Successfully built spiceypy
Installing collected packages: spiceypy
Successfully installed spiceypy-4.0.0
In [3]:
"ASTROPY"
!pip install astroquery
from astroquery.jplhorizons import Horizons

from astropy.time import Time
from astropy.table import Table
from astropy.coordinates import SkyCoord, EarthLocation, AltAz, get_sun
import astropy.units as u
Collecting astroquery
  Downloading https://files.pythonhosted.org/packages/1b/f8/4690523783691ed816b3469c3ec611af3798594d37ade510dd918d59f57e/astroquery-0.4.1.tar.gz (6.5MB)
     |████████████████████████████████| 6.5MB 4.8MB/s 
Requirement already satisfied: numpy in /usr/local/lib/python3.6/dist-packages (from astroquery) (1.19.5)
Requirement already satisfied: astropy>=3.1 in /usr/local/lib/python3.6/dist-packages (from astroquery) (4.1)
Requirement already satisfied: requests>=2.4.3 in /usr/local/lib/python3.6/dist-packages (from astroquery) (2.23.0)
Collecting keyring>=4.0
  Downloading https://files.pythonhosted.org/packages/d0/a0/20e656cd1e2313af619e382782bd47b5f77a3f33d81992554f3aac56e90d/keyring-21.8.0-py3-none-any.whl
Requirement already satisfied: beautifulsoup4>=4.3.2 in /usr/local/lib/python3.6/dist-packages (from astroquery) (4.6.3)
Requirement already satisfied: html5lib>=0.999 in /usr/local/lib/python3.6/dist-packages (from astroquery) (1.0.1)
Requirement already satisfied: six in /usr/local/lib/python3.6/dist-packages (from astroquery) (1.15.0)
Requirement already satisfied: urllib3!=1.25.0,!=1.25.1,<1.26,>=1.21.1 in /usr/local/lib/python3.6/dist-packages (from requests>=2.4.3->astroquery) (1.24.3)
Requirement already satisfied: certifi>=2017.4.17 in /usr/local/lib/python3.6/dist-packages (from requests>=2.4.3->astroquery) (2020.12.5)
Requirement already satisfied: idna<3,>=2.5 in /usr/local/lib/python3.6/dist-packages (from requests>=2.4.3->astroquery) (2.10)
Requirement already satisfied: chardet<4,>=3.0.2 in /usr/local/lib/python3.6/dist-packages (from requests>=2.4.3->astroquery) (3.0.4)
Requirement already satisfied: importlib-metadata>=1; python_version < "3.8" in /usr/local/lib/python3.6/dist-packages (from keyring>=4.0->astroquery) (3.3.0)
Collecting jeepney>=0.4.2; sys_platform == "linux"
  Downloading https://files.pythonhosted.org/packages/51/b0/a6ea72741aaac3f37fb96d195e4ee576a103c4c04e279bc6b446a70960e1/jeepney-0.6.0-py3-none-any.whl (45kB)
     |████████████████████████████████| 51kB 5.7MB/s 
Collecting SecretStorage>=3.2; sys_platform == "linux"
  Downloading https://files.pythonhosted.org/packages/63/a2/a6d9099b14eb5dbbb04fb722d2b5322688f8f99b471bdf2097e33efa8091/SecretStorage-3.3.0-py3-none-any.whl
Requirement already satisfied: webencodings in /usr/local/lib/python3.6/dist-packages (from html5lib>=0.999->astroquery) (0.5.1)
Requirement already satisfied: zipp>=0.5 in /usr/local/lib/python3.6/dist-packages (from importlib-metadata>=1; python_version < "3.8"->keyring>=4.0->astroquery) (3.4.0)
Requirement already satisfied: typing-extensions>=3.6.4; python_version < "3.8" in /usr/local/lib/python3.6/dist-packages (from importlib-metadata>=1; python_version < "3.8"->keyring>=4.0->astroquery) (3.7.4.3)
Collecting cryptography>=2.0
  Downloading https://files.pythonhosted.org/packages/c9/de/7054df0620b5411ba45480f0261e1fb66a53f3db31b28e3aa52c026e72d9/cryptography-3.3.1-cp36-abi3-manylinux2010_x86_64.whl (2.6MB)
     |████████████████████████████████| 2.6MB 41.1MB/s 
Requirement already satisfied: cffi>=1.12 in /usr/local/lib/python3.6/dist-packages (from cryptography>=2.0->SecretStorage>=3.2; sys_platform == "linux"->keyring>=4.0->astroquery) (1.14.4)
Requirement already satisfied: pycparser in /usr/local/lib/python3.6/dist-packages (from cffi>=1.12->cryptography>=2.0->SecretStorage>=3.2; sys_platform == "linux"->keyring>=4.0->astroquery) (2.20)
Building wheels for collected packages: astroquery
  Building wheel for astroquery (setup.py) ... done
  Created wheel for astroquery: filename=astroquery-0.4.1-cp36-none-any.whl size=3831873 sha256=e2d7f0574202bd140a3313f4285985432affa63685448eb04f69344902df7130
  Stored in directory: /root/.cache/pip/wheels/88/f8/b7/a254cd96e808f708bc0b7d755a8e095c56fbbe94099d7b464f
Successfully built astroquery
Installing collected packages: jeepney, cryptography, SecretStorage, keyring, astroquery
Successfully installed SecretStorage-3.3.0 astroquery-0.4.1 cryptography-3.3.1 jeepney-0.6.0 keyring-21.8.0
In [4]:
"MAPITAS"
!apt-get install libgeos-3.5.0
!apt-get install libgeos-dev
!pip install https://github.com/matplotlib/basemap/archive/master.zip
!pip install pyproj==1.9.6
from mpl_toolkits.basemap import Basemap
Reading package lists... Done
Building dependency tree       
Reading state information... Done
E: Unable to locate package libgeos-3.5.0
E: Couldn't find any package by glob 'libgeos-3.5.0'
E: Couldn't find any package by regex 'libgeos-3.5.0'
Reading package lists... Done
Building dependency tree       
Reading state information... Done
Suggested packages:
  libgdal-doc
The following NEW packages will be installed:
  libgeos-dev
0 upgraded, 1 newly installed, 0 to remove and 16 not upgraded.
Need to get 73.1 kB of archives.
After this operation, 486 kB of additional disk space will be used.
Get:1 http://archive.ubuntu.com/ubuntu bionic/universe amd64 libgeos-dev amd64 3.6.2-1build2 [73.1 kB]
Fetched 73.1 kB in 0s (187 kB/s)
Selecting previously unselected package libgeos-dev.
(Reading database ... 145483 files and directories currently installed.)
Preparing to unpack .../libgeos-dev_3.6.2-1build2_amd64.deb ...
Unpacking libgeos-dev (3.6.2-1build2) ...
Setting up libgeos-dev (3.6.2-1build2) ...
Processing triggers for man-db (2.8.3-2ubuntu0.1) ...
Collecting https://github.com/matplotlib/basemap/archive/master.zip
  Downloading https://github.com/matplotlib/basemap/archive/master.zip (133.1MB)
     |████████████████████████████████| 133.1MB 85kB/s 
Requirement already satisfied: matplotlib!=3.0.1,>=1.0.0 in /usr/local/lib/python3.6/dist-packages (from basemap==1.2.2+dev) (3.2.2)
Requirement already satisfied: numpy>=1.2.1 in /usr/local/lib/python3.6/dist-packages (from basemap==1.2.2+dev) (1.19.5)
Collecting pyproj>=1.9.3
  Downloading https://files.pythonhosted.org/packages/e4/ab/280e80a67cfc109d15428c0ec56391fc03a65857b7727cf4e6e6f99a4204/pyproj-3.0.0.post1-cp36-cp36m-manylinux2010_x86_64.whl (6.4MB)
     |████████████████████████████████| 6.5MB 5.1MB/s 
Collecting pyshp>=1.2.0
  Downloading https://files.pythonhosted.org/packages/38/85/fbf87e7aa55103e0d06af756bdbc15cf821fa580414c23142d60a35d4f85/pyshp-2.1.3.tar.gz (219kB)
     |████████████████████████████████| 225kB 48.5MB/s 
Requirement already satisfied: six in /usr/local/lib/python3.6/dist-packages (from basemap==1.2.2+dev) (1.15.0)
Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib!=3.0.1,>=1.0.0->basemap==1.2.2+dev) (2.4.7)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.6/dist-packages (from matplotlib!=3.0.1,>=1.0.0->basemap==1.2.2+dev) (0.10.0)
Requirement already satisfied: python-dateutil>=2.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib!=3.0.1,>=1.0.0->basemap==1.2.2+dev) (2.8.1)
Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib!=3.0.1,>=1.0.0->basemap==1.2.2+dev) (1.3.1)
Requirement already satisfied: certifi in /usr/local/lib/python3.6/dist-packages (from pyproj>=1.9.3->basemap==1.2.2+dev) (2020.12.5)
Building wheels for collected packages: basemap, pyshp
  Building wheel for basemap (setup.py) ... done
  Created wheel for basemap: filename=basemap-1.2.2+dev-cp36-cp36m-linux_x86_64.whl size=121757244 sha256=440f16a28379a65e30dcec1f785e0b2a1e06e98d179c0d533da6a2e2b9512c68
  Stored in directory: /tmp/pip-ephem-wheel-cache-sex59hoy/wheels/98/4a/fc/ce719b75d97e646645c225f3332b1b217536100314922e9572
  Building wheel for pyshp (setup.py) ... done
  Created wheel for pyshp: filename=pyshp-2.1.3-cp36-none-any.whl size=37264 sha256=db1f58dab753afd0b7ba74d1c594814207b28ca2734fbd71d723ded0fa73c198
  Stored in directory: /root/.cache/pip/wheels/76/2b/d4/53e6b9a0fb0a9f9f29664cf82605af8bc81d5ab44d987896dd
Successfully built basemap pyshp
Installing collected packages: pyproj, pyshp, basemap
Successfully installed basemap-1.2.2+dev pyproj-3.0.0.post1 pyshp-2.1.3
Collecting pyproj==1.9.6
  Downloading https://files.pythonhosted.org/packages/26/8c/1da0580f334718e04f8bbf74f0515a7fb8185ff96b2560ce080c11aa145b/pyproj-1.9.6.tar.gz (2.8MB)
     |████████████████████████████████| 2.8MB 4.3MB/s 
Building wheels for collected packages: pyproj
  Building wheel for pyproj (setup.py) ... done
  Created wheel for pyproj: filename=pyproj-1.9.6-cp36-cp36m-linux_x86_64.whl size=3702125 sha256=e8fe3a8c565121753dd7404b22b3f2783d94599be23f9f66d07b91e1f91e73b7
  Stored in directory: /root/.cache/pip/wheels/02/cd/b1/a2d6430f74c7a778a43d62f78bec109ca69c732dc9b929142a
Successfully built pyproj
Installing collected packages: pyproj
  Found existing installation: pyproj 3.0.0.post1
    Uninstalling pyproj-3.0.0.post1:
      Successfully uninstalled pyproj-3.0.0.post1
Successfully installed pyproj-1.9.6
In [5]:
"TQDM"
!pip install tqdm
from tqdm import tqdm
Requirement already satisfied: tqdm in /usr/local/lib/python3.6/dist-packages (4.41.1)

Efemérides

In [6]:
"UNIDAD ASTRONÓMICA"
AU = (1*u.au).to(u.m).value

"TIEMPOS"
# Se definen los tiempos desde que está a 12Rt, pasando por su máximo
# acercamiento en 6Rt y volviendo a 12Rt
t_ini = Time('2029-04-13 19:00:00',format='iso')
t_max = Time('2029-04-13 22:05:00',format='iso') #Máximo Acercamiento
# t_end = Time('2029-04-13 23:00:00',format='iso') #Final. Step de 5min
t_end = Time('2029-04-14 06:00:00',format='iso')

timesyy = {'start':t_ini.value, 'stop':t_end.value, 'step':'20m'}


print("(Ignorar advertencia de ERFA si sale, es porque está muy a futuro)")
(Ignorar advertencia de ERFA si sale, es porque está muy a futuro)
WARNING: ErfaWarning: ERFA function "dtf2d" yielded 1 of "dubious year (Note 6)" [astropy._erfa.core]
WARNING: ErfaWarning: ERFA function "d2dtf" yielded 1 of "dubious year (Note 5)" [astropy._erfa.core]
In [7]:
"EFEMÉRIDES"
#Efemérides de Apofis entre los tiempos estipulados
apophis = Horizons(id='99942',epochs=timesyy) 
ephemeris = apophis.ephemerides()
vectors = apophis.vectors()

#Coordenadas Ecuatoriales del Asteroide
ECUcoords = SkyCoord(ephemeris['RA'],ephemeris['DEC']
                     ,ephemeris['delta'],frame='gcrs')

Nt = np.shape(ephemeris)[0]
print(Nt)
34
In [8]:
# ephemeris.to_pandas()
plt.plot(ephemeris.to_pandas()['datetime_jd']-2.462240e+06,ephemeris.to_pandas()['delta']*23500)
Out[8]:
[<matplotlib.lines.Line2D at 0x7ff8d8138fd0>]

Coloreando el Mapita

In [9]:
def get_h(T,Lon,Lat,RA,Dec,r):
  """
  Teniendo las coordenadas ecuatoriales y la ubicación de un observador,
  se saca el ángulo de elevación.
  """
  G = T.sidereal_time('apparent', 'greenwich').to(u.deg)
  LonR = Lon*u.deg + G

  #Vector Posición Objeto
  Rx = r*np.cos(Dec*u.deg)*np.cos(RA*u.deg)
  Ry = r*np.cos(Dec*u.deg)*np.sin(RA*u.deg)
  Rz = r*np.sin(Dec*u.deg)
  R = np.array([Rx,Ry,Rz])

  #Vector Posición Observador
  RLx = np.cos(Lat*u.deg)*np.cos(LonR)
  RLy = np.cos(Lat*u.deg)*np.sin(LonR)
  RLz = np.sin(Lat*u.deg)
  RL = np.array([RLx,RLy,RLz])

  #Resta de Vectores y Producto Escalar Inverso
  R2 = R-RL
  H = np.rad2deg(np.arccos(np.dot(RL,R2)/(np.linalg.norm(RL)*np.linalg.norm(R2))))
  H = 90 - H
  return H
In [10]:
"SACO LA ALTURA DEL ASTEROIDE EN TODAS LAS UBICACIONES (SE DEMORA COMO 2 MINUTOS)"
nlon = 64
nlat = 36
Lons = np.linspace(-180,180,nlon)
Lats = np.linspace(-90,90,nlat)

#Creo el Grid de Latitudes, Longitudes y fechas
Alts = np.zeros((nlon,nlat,Nt))

#Este siempre se demora ratico
for t in tqdm(range(Nt)):
  TIME = Time(ephemeris['datetime_jd'][t],format='jd')
  RA = ephemeris['RA'][t]
  DEC = ephemeris['DEC'][t]
  DELTA = ephemeris['delta'][t]*u.au.to(u.earthRad)
  for i in range(nlon):
    for j in range(nlat):
      LON = Lons[i]
      LAT = Lats[j]
      Alts[i,j,t] = get_h(TIME,LON,LAT,RA,DEC,DELTA)
  0%|          | 0/34 [00:00<?, ?it/s]WARNING: ErfaWarning: ERFA function "utcut1" yielded 1 of "dubious year (Note 3)" [astropy._erfa.core]
WARNING: ErfaWarning: ERFA function "utctai" yielded 1 of "dubious year (Note 3)" [astropy._erfa.core]
100%|██████████| 34/34 [01:38<00:00,  2.90s/it]
In [11]:
plt.figure(figsize=(10,6))
plt.contourf(Lons,Lats,Alts[:,:,0].T,levels=np.linspace(-90,90,100),cmap='inferno')
plt.colorbar()
plt.title("Alturas de Apophis en "+str(ephemeris['datetime_str'][0]))
Out[11]:
Text(0.5, 1.0, 'Alturas de Apophis en 2029-Apr-13 19:00')
In [12]:
def filterNeg(arr):
  """
  Recibe un array en 3D y le quita los valores negativos
  """
  a,b,c = np.shape(arr)
  arr1 = np.reshape(arr,a*b*c)
  arr1[np.where(arr1 < 0)] = 0

  return np.reshape(arr1,(a,b,c))/90
In [13]:
AltsPos = filterNeg(Alts)

plt.figure(figsize=(10,6))
plt.contourf(Lons,Lats,AltsPos[:,:,0].T,levels=np.linspace(0,1,100),cmap='inferno')
plt.colorbar()
plt.title("Alturas Positivas de Apophis en "+str(ephemeris['datetime_str'][0]))
Out[13]:
Text(0.5, 1.0, 'Alturas Positivas de Apophis en 2029-Apr-13 19:00')
In [14]:
"ALTURAS DEL SOL (TAMBIÉN SE DEMORA COMO 2 MINUTOS)"
AltsSun = np.zeros((nlon,nlat,Nt))
FilterSun = np.zeros((nlon,nlat,Nt))

#Defino una función de caída para filtrar con la altura del Sol
def sigmoid(x):
  return 1 / (1 + np.exp(0.5*(x+10)))

for t in tqdm(range(Nt)):
  TIME = Time(ephemeris['datetime_jd'][t],format='jd')
  SUN = get_sun(TIME)
  Sun_RA = SUN.ra.value
  Sun_DEC = SUN.dec.value
  Sun_DIST = 1*u.au.to(u.earthRad)
  for i in range(nlon):
    for j in range(nlat):
      LON = Lons[i]
      LAT = Lats[j]

      Sun_H = get_h(TIME,LON,LAT,
                             Sun_RA,Sun_DEC,Sun_DIST)
      AltsSun[i,j,t] = Sun_H
      FilterSun[i,j,t] = sigmoid(Sun_H)
  0%|          | 0/34 [00:00<?, ?it/s]WARNING: ErfaWarning: ERFA function "utctai" yielded 1 of "dubious year (Note 3)" [astropy._erfa.core]
WARNING: ErfaWarning: ERFA function "taiutc" yielded 1 of "dubious year (Note 4)" [astropy._erfa.core]
WARNING: ErfaWarning: ERFA function "utcut1" yielded 1 of "dubious year (Note 3)" [astropy._erfa.core]
100%|██████████| 34/34 [01:39<00:00,  2.92s/it]
In [15]:
plt.figure(figsize=(10,6))
plt.contourf(Lons,Lats,AltsSun[:,:,0].T,levels=np.linspace(-90,90,100),cmap='inferno')
plt.colorbar()
plt.title("Alturas del Sol en "+str(ephemeris['datetime_str'][0]))
Out[15]:
Text(0.5, 1.0, 'Alturas del Sol en 2029-Apr-13 19:00')
In [16]:
AltsPosSun = AltsPos * FilterSun

plt.figure(figsize=(10,6))
plt.contourf(Lons,Lats,AltsPosSun[:,:,0].T,levels=np.linspace(0,1,100),cmap='inferno')
plt.colorbar()
plt.title("Alturas de Apophis donde el Sol no Molesta \n en "+str(ephemeris['datetime_str'][0]))
plt.show()

Plot

In [17]:
ephemeris.to_pandas().columns
Out[17]:
Index(['targetname', 'datetime_str', 'datetime_jd', 'H', 'G', 'solar_presence',
       'flags', 'RA', 'DEC', 'RA_app', 'DEC_app', 'RA_rate', 'DEC_rate', 'AZ',
       'EL', 'AZ_rate', 'EL_rate', 'sat_X', 'sat_Y', 'sat_PANG',
       'siderealtime', 'airmass', 'magextinct', 'V', 'surfbright',
       'illumination', 'illum_defect', 'sat_sep', 'sat_vis', 'ang_width',
       'PDObsLon', 'PDObsLat', 'PDSunLon', 'PDSunLat', 'SubSol_ang',
       'SubSol_dist', 'NPole_ang', 'NPole_dist', 'EclLon', 'EclLat', 'r',
       'r_rate', 'delta', 'delta_rate', 'lighttime', 'vel_sun', 'vel_obs',
       'elong', 'elongFlag', 'alpha', 'lunar_elong', 'lunar_illum',
       'sat_alpha', 'sunTargetPA', 'velocityPA', 'OrbPlaneAng',
       'constellation', 'TDB-UT', 'ObsEclLon', 'ObsEclLat', 'NPole_RA',
       'NPole_DEC', 'GlxLon', 'GlxLat', 'solartime', 'earth_lighttime',
       'RA_3sigma', 'DEC_3sigma', 'SMAA_3sigma', 'SMIA_3sigma', 'Theta_3sigma',
       'Area_3sigma', 'RSS_3sigma', 'r_3sigma', 'r_rate_3sigma',
       'SBand_3sigma', 'XBand_3sigma', 'DoppDelay_3sigma', 'true_anom',
       'hour_angle', 'alpha_true', 'PABLon', 'PABLat'],
      dtype='object')
In [18]:
def get_lon(TIME,RA):
  """
  Saca la Longitud que estará bajo el asteroide en cada momento
  """
  T = Time(TIME,format='jd')
  G = T.sidereal_time('apparent', 'greenwich').to(u.deg)
  return RA - G.value

maxlons = get_lon(ephemeris['datetime_jd'],ephemeris['RA'])
maxlats = ephemeris['DEC']

for i in range(len(maxlons)):
  if maxlons[i] < -180: maxlons[i] += 360

# for i in range(1,len(maxlons)):
#   if maxlons[i] > maxlons[i-1]:

i = -15
maxlons[i] = np.nan
maxlons = np.array(np.insert(maxlons,i+1,[179]))
maxlats = np.array(np.insert(maxlats,i+1,[maxlats[i-1]]))
WARNING: ErfaWarning: ERFA function "utcut1" yielded 34 of "dubious year (Note 3)" [astropy._erfa.core]
WARNING: ErfaWarning: ERFA function "utctai" yielded 34 of "dubious year (Note 3)" [astropy._erfa.core]
In [19]:
maxlons[-20:],maxlats[-20:]
Out[19]:
(array([-151.82111506, -161.63334432, -170.52351359, -178.71490285,
                  nan,  179.        ,  166.38950894,  159.47785967,
         152.8284904 ,  146.39158112,  140.12887184,  134.01059256,
         128.01329328,  122.118244  ,  116.31038472,  110.57742577,
         104.90927648,   99.29754719,   93.7352679 ,   88.2165386 ]),
 array([36.23366, 34.96099, 33.76112, 32.65895, 31.65889, 32.65895,
        30.75613, 29.94222, 29.20775, 28.54357, 27.94132, 27.39355,
        26.89374, 26.43627, 26.01623, 25.62941, 25.27216, 24.94132,
        24.63413, 24.3482 ]))
In [20]:
ephemeris['datetime_str'][-12] #Paso Máximo AQUI
Out[20]:
'2029-Apr-14 02:20'
In [21]:
from mpl_toolkits.axes_grid1.inset_locator import inset_axes

t = -1
fig=plt.figure(figsize=(10,8))
ax = plt.axes()
ax.clear()
plot = ax.contourf(Lons,Lats,AltsPosSun[:,:,t].T*90,levels=np.linspace(-0.1,1,100)*90,cmap='Blues')
m = Basemap(ax=ax,projection='cyl', lon_0 = 0, lat_0 = 0)
m.drawcoastlines()
m.drawparallels(np.arange(-90,90,30),labels=[1,1,0,1])
m.drawmeridians(np.arange(-180,180,30),labels=[1,1,0,1])
ax.set_title("Trayectoria y Elevaciones para el Asteroide 99942 Apophis \n"+str(ephemeris['datetime_str'][t]))


ax.plot(maxlons,maxlats,'w-',lw=5)
ax.plot(maxlons[t],maxlats[t],'wo',ms=12)

ax.plot(maxlons,maxlats,'b-',lw=2)
ax.plot(maxlons[t],maxlats[t],'bo',ms=9)


axins = inset_axes(ax,
                   width="5%",  # width = 5% of parent_bbox width
                   height="100%",  # height : 50%
                   loc='lower left',
                   bbox_to_anchor=(1.07, 0., 1, 1),
                   bbox_transform=ax.transAxes,
                   borderpad=0,
                   )

cbar = fig.colorbar(plot, cax=axins)
cbar.ax.get_yaxis().labelpad = 15
cbar.ax.set_ylabel('Elevación sobre el Horizonte [°]', rotation=270)
plt.show()
In [22]:
def frame(t):
    ax.clear()
    plot = ax.contourf(Lons,Lats,AltsPosSun[:,:,t].T*90,levels=np.linspace(-0.1,1,100)*90,cmap='Blues')
    m = Basemap(ax=ax,projection='cyl', lon_0 = 0, lat_0 = 0)
    m.drawcoastlines()
    m.drawparallels(np.arange(-90,90,30),labels=[1,1,0,1])
    m.drawmeridians(np.arange(-180,180,30),labels=[1,1,0,1])
    ax.set_title("Mapa de Contorno de Elevaciones y Trayectoria para el Asteroide 99942 Apophis \n"+str(ephemeris['datetime_str'][t]))

    ax.plot(maxlons,maxlats,'w-',lw=5)
    ax.plot(maxlons[t],maxlats[t],'wo',ms=12)

    ax.plot(maxlons,maxlats,'b-',lw=2)
    ax.plot(maxlons[t],maxlats[t],'bo',ms=9)


    axins = inset_axes(ax,
                      width="5%",  # width = 5% of parent_bbox width
                      height="100%",  # height : 50%
                      loc='lower left',
                      bbox_to_anchor=(1.07, 0., 1, 1),
                      bbox_transform=ax.transAxes,
                      borderpad=0,
                      )

    cbar = fig.colorbar(plot, cax=axins)
    cbar.ax.get_yaxis().labelpad = 15
    cbar.ax.set_ylabel('Elevación sobre el Horizonte [°]', rotation=270)
    return plot

Animación

In [23]:
fig=plt.figure(figsize=(10,8))
ax = plt.axes()

anim = animation.FuncAnimation(fig, frame, frames=Nt, blit=False, repeat=True)
In [24]:
Writer = animation.writers['ffmpeg']
writer = Writer(fps=5, metadata=dict(artist='Agustín-Vallejo'), bitrate=1800)
anim.save("Apophis3.mp4",writer=writer)

from google.colab import files
files.download("Apophis3.mp4")
In [25]:
anim
Out[25]:
In [26]:
from google.colab import drive
drive.mount('/content/drive')
Mounted at /content/drive

FAST

In [ ]:
"UNIDAD ASTRONÓMICA"
AU = (1*u.au).to(u.m).value

"TIEMPOS"
# Se definen los tiempos desde que está a 12Rt, pasando por su máximo
# acercamiento en 6Rt y volviendo a 12Rt
t_ini = Time('2029-04-13 19:00:00',format='iso')
t_max = Time('2029-04-13 22:05:00',format='iso') #Máximo Acercamiento
t_end = Time('2029-04-14 12:00:00',format='iso')
timesyy = {'start':t_ini.value, 'stop':t_end.value, 'step':'15m'}


print("(Ignorar advertencia de ERFA si sale, es porque está muy a futuro)")
(Ignorar advertencia de ERFA si sale, es porque está muy a futuro)
WARNING: ErfaWarning: ERFA function "dtf2d" yielded 1 of "dubious year (Note 6)" [astropy._erfa.core]
WARNING: ErfaWarning: ERFA function "d2dtf" yielded 1 of "dubious year (Note 5)" [astropy._erfa.core]
In [ ]:
"EFEMÉRIDES"
#Efemérides de Apofis entre los tiempos estipulados
arecibo = { 'lon': 293.24692,
            'lat': 0.949577,
            'elevation': +0.312734}

fast   = { 'lon': 106.85667,
            'lat': 25.653055,
            'elevation': +0.312734}
apophis = Horizons(id='99942',location=fast,epochs=timesyy) 
ephemeris = apophis.ephemerides()
apophis.vectors()
In [ ]:
ephemeris.to_pandas().columns
Out[ ]:
Index(['targetname', 'datetime_str', 'datetime_jd', 'H', 'G', 'solar_presence',
       'flags', 'RA', 'DEC', 'RA_app', 'DEC_app', 'RA_rate', 'DEC_rate', 'AZ',
       'EL', 'AZ_rate', 'EL_rate', 'sat_X', 'sat_Y', 'sat_PANG',
       'siderealtime', 'airmass', 'magextinct', 'V', 'surfbright',
       'illumination', 'illum_defect', 'sat_sep', 'sat_vis', 'ang_width',
       'PDObsLon', 'PDObsLat', 'PDSunLon', 'PDSunLat', 'SubSol_ang',
       'SubSol_dist', 'NPole_ang', 'NPole_dist', 'EclLon', 'EclLat', 'r',
       'r_rate', 'delta', 'delta_rate', 'lighttime', 'vel_sun', 'vel_obs',
       'elong', 'elongFlag', 'alpha', 'lunar_elong', 'lunar_illum',
       'sat_alpha', 'sunTargetPA', 'velocityPA', 'OrbPlaneAng',
       'constellation', 'TDB-UT', 'ObsEclLon', 'ObsEclLat', 'NPole_RA',
       'NPole_DEC', 'GlxLon', 'GlxLat', 'solartime', 'earth_lighttime',
       'RA_3sigma', 'DEC_3sigma', 'SMAA_3sigma', 'SMIA_3sigma', 'Theta_3sigma',
       'Area_3sigma', 'RSS_3sigma', 'r_3sigma', 'r_rate_3sigma',
       'SBand_3sigma', 'XBand_3sigma', 'DoppDelay_3sigma', 'true_anom',
       'hour_angle', 'alpha_true', 'PABLon', 'PABLat'],
      dtype='object')
In [ ]:
ephemeris['horas'] = 24*(ephemeris['datetime_jd']-ephemeris['datetime_jd'][0])
In [ ]:
H1 = ephemeris['horas'][np.where(ephemeris['EL'] >= 50)][0]
H2 = ephemeris['horas'][np.where(ephemeris['EL'] >= 50)][-1]
print(H1,H2)
7.749999985098839 12.49999999254942
In [ ]:
ephemeris[['datetime_str','RA','DEC','EL']][np.where(ephemeris['EL'] >= 49)].to_pandas()
Out[ ]:
datetime_str RA DEC EL
0 2029-Apr-14 02:30 31.93984 28.66451 49.708346
1 2029-Apr-14 02:45 30.61484 28.25000 54.108386
2 2029-Apr-14 03:00 29.40116 27.85451 58.447295
3 2029-Apr-14 03:15 28.28617 27.47722 62.730263
4 2029-Apr-14 03:30 27.25914 27.11721 66.961798
5 2029-Apr-14 03:45 26.31091 26.77352 71.145817
6 2029-Apr-14 04:00 25.43359 26.44523 75.285713
7 2029-Apr-14 04:15 24.62040 26.13143 79.384396
8 2029-Apr-14 04:30 23.86545 25.83126 83.444257
9 2029-Apr-14 04:45 23.16361 25.54391 87.466149
10 2029-Apr-14 05:00 22.51039 25.26864 88.534239
11 2029-Apr-14 05:15 21.90182 25.00474 84.581698
12 2029-Apr-14 05:30 21.33441 24.75159 80.656515
13 2029-Apr-14 05:45 20.80505 24.50859 76.760875
14 2029-Apr-14 06:00 20.31096 24.27519 72.894006
15 2029-Apr-14 06:15 19.84964 24.05091 69.055091
16 2029-Apr-14 06:30 19.41885 23.83529 65.243420
17 2029-Apr-14 06:45 19.01656 23.62790 61.458418
18 2029-Apr-14 07:00 18.64092 23.42837 57.699636
19 2029-Apr-14 07:15 18.29023 23.23634 53.966760
20 2029-Apr-14 07:30 17.96295 23.05148 50.259602
In [ ]:
fig,ax = plt.subplots(figsize=(8,5))

ax.set_title("Efemérides de 99942 Apophis desde FAST (25.6°N - 106.8°E)")
ax.plot(ephemeris['horas'],ephemeris['EL'],label='Elevación [°]',lw=3,color='b')
# plt.plot(ephemeris.to_pandas()['horas'],ephemeris.to_pandas()['DEC'],label='Declinación [°]',lw=3)
# plt.plot(ephemeris.to_pandas()['horas'],ephemeris.to_pandas()['elong'],label='Elongación [°]',lw=3)

ax.fill_between([H1-0.2,H2],0,100,alpha=0.2,color='b')
ax.text(8.5,45,"Apophis Visible \n desde FAST",color='b')
# ax.hlines(50,-10,50,'b',linestyles='--')

ax.set_ylabel("Elevación del Apophis (°)",color='b',size=14)

plt.xlabel("UT el 2029-04-14",size=15)
plt.xticks([0,5,10,15],['-5:00','0:00','5:00','10:00'],size=15)
plt.yticks(size=12)
plt.ylim([0,90])
plt.xlim([0,ephemeris['horas'][-1]])
ax.grid()

ax2=ax.twinx()
ax2.plot(ephemeris['horas'],ephemeris['delta']*23500,label='Distancia [Rt]',lw=3,color='r')
ax.plot(ephemeris['horas'],ephemeris['delta']*0-1,label='Distancia [Rt]',lw=3,color='r')

ax2.set_ylabel("Distancia al Apophis ($R_t$)",color='r',size=14)


ax.legend(loc=2)
plt.show()